rm(list=ls(all=TRUE))
load("~/Post-doc/Data/Total Merged Data File (July 24 2023).RData")
str(CC.TotalData)
Data for mean velocity, angular velocity, turn and heading angles each will be analyzed by a mixed model ANCOVA with chemical concentration, and light as fixed (treatment) effects, flow velocity as the co-variate and replicate as a random (blocking affect) with individual paths as the unit of analysis.
##To run ANCOVA in R load the following packages:
library(car)
library(compute.es)
library(effects)
library(ggplot2)
library(multcomp)
library(pastecs)
library(WRS2)
##If you are using lm or aov make sure that you set the contrasts using the "contrasts" function before doing either aov or lm. R uses non-orthogonal contrasts by default which can mess everything up in an ANCOVA. If you want to set orthogonal contrasts use:
contrasts(dataname$factorvariable)=contr.poly()# of levels, i.e. 3) ##example
contrasts(CC.TotalData$Chlorophyll)=contr.poly(7)
##then run your model as
model.1=aov(dv~covariate+factorvariable, data=dataname) ## example
mod1=aov(vel.flow~Flow.rate+Chlorophyll, data=CC.TotalData)
##To view the model use:
Anova(model.1, type="III") ## example
Anova(mod1, type="III")
##Make sure you use capital "A" Anova here and not anova. This will give results using type III SS.
summary.lm(model.1) ##will give another summary and includes the R-sq. output. Example
summary.lm(mod1)
posth=glht(model.1, linfct=mcp(factorvariable="Tukey")) ##gives the post-hoc Tukey analysis. Example
posth=glht(mod1, linfct=mcp(Chlorophyll="Tukey"))
summary(posth) ##shows the output in a nice format.
##If you want to test for homogeneity of regression slopes you can also include an interaction term for the IV and covariate. That would be:
model=aov(dv~covariate+IV+covariate:IV, data=dataname) ## example
mod2=aov(vel.flow~Flow.rate+Chlorophyll+Flow.rate:Chlorophyll, data=CC.TotalData) ## not right?
##If the interaction term is significant then you do not have homogeneity.
Stats on each level
##Looking at other stats we could run
library(nlme)
library(lme4)
library(effects)
## example
d <- data.frame(
Y = rnorm(48), ### e.g. velocity
subject = factor(rep(1:12, 4)), ### e.g. D_V_T
A = factor(rep(1:2, each=24)), #### e.g. Light
B = factor(rep(rep(1:2, each=12), 2))) ### e.g. Flow.rate
summary(aov(Y ~ A*B + Error(subject/(A*B)), data=d)) # Standard repeated measures ANOVA
mode1 <- lme(Y~A*B, random = ~1|subject, data=d) ## nlme
lmer(Y ~ A*B + (1|subject) + (1|A:subject) + (1|B:subject), data=d) ### lme4
######### our data using nlme package #########
mode1 <- lme(pitch.perfect~Chlorophyll*Light, random= ~1|D_V_T, data=CC.TotalData)
mode2 <- lme(pitch.perfect~Light, random= ~1|D_V_T, data=CC.TotalData)
anova(mode1, mode2) #if it shows only minor improvement, no need to include missing factor
anova(mode2)
######### our data using l4me package #########
model1 <- lmer(vel.flow ~ Flow.rate*Chlorophyll + (1|D_V_T) + (1|Flow.rate:D_V_T) + (1|Chlorophyll:D_V_T), data=CC.TotalData) ###
model2 <- lmer(vel.flow ~ Flow.rate + (1|D_V_T) + (1|Flow.rate:D_V_T), data=CC.TotalData)
anova(model1, model2) #if it shows only minor improvement, no need to include missing factor
anova(model1)
######### plotting effect of the models ##########
plot(effect("turn.angle*Flow.rate",se=TRUE, confidence.level=.95, mode2))
###### interaction plot with Lattice extra #########
library(latticeExtra)
xyplot(CC.TotalData$vel.flow~CC.TotalData$turn.angle|CC.TotalData$Light*CC.TotalData$Flow.rate, main="Scatterplots by Light and Flow Rate",
ylab="Vel.flow", xlab="Turn Angle") +
layer(panel.ablineq(lm(y ~ x), r.sq = TRUE,
at = 0.75, adj=1:5), style = 4)
############## interaction plots ###########################
interaction.plot(CC.TotalData$vel.flow, ## x axis
CC.TotalData$Flow.rate, ## factor for lines
CC.TotalData$turn.angle, ## y axis
fun = mean,
type = c("b"), legend = TRUE,
trace.label = "Flow Rate",
ylim = c(-1,180),
main = "Lights On",
xlab = "vel.flow",
ylab = "turn angle",
col = CC.TotalData$Flow.rate, pch = c(1:9),
axes = TRUE)
#library(tidyr)
#library(dplyr)
#CC.TotalData %>% mutate_all(~replace(., is.nan(.), 0))
CC.TotalData[is.nan(CC.TotalData)] <- 0
range(CC.TotalData$log.vel.flow)
##########################################################
rm(mode1)
rm(mode2)
mod1 <- lm(vel.flow~Flow.rate*Chlorophyll, data = CC.TotalData)
summary(mode2)
Kernel Density Plots
Turn Angles first…
################################ Chlorophyll ########################
levels(CC.TotalData$Chlorophyll)
[1] "No Chlorophyll" "Low Chlorophyll" "Medium Chlorophyll"
[4] "High Chlorophyll"
Now looking at velocity in relation to flow instead….
### Kernel density plots ################### Swimming velocity
##range(CC.TotalData$vel.flow)
## plot(log.vel.flow)
##log.vel.flow <- log10(CC.TotalData$vel.flow[CC.TotalData$vel.flow > 0])
#str(CC.TotalData)
# vels <- (CC.TotalData$v[CC.TotalData$D_V_T==ind[i]]) ## vels <- log10(vels[vels>0]) ## from notebook01, could be useful to try this instead
################################ Flow Rate ########################
levels(CC.TotalData$Flow.rate)
[1] "0" "0.6" "3" "5.9" "8.9"
Now for headings…..
rm(kd1, kd2, kd3, kd4, kd5, kd6, kd7, kd8)
Warning in rm(kd1, kd2, kd3, kd4, kd5, kd6, kd7, kd8) :
object 'kd6' not found
Warning in rm(kd1, kd2, kd3, kd4, kd5, kd6, kd7, kd8) :
object 'kd7' not found
Warning in rm(kd1, kd2, kd3, kd4, kd5, kd6, kd7, kd8) :
object 'kd8' not found
Now for pitch….